Solution of the Percus-Yevick equation for hard hyperspheres in even dimensions 
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We solve the Percus-Yevick equation in even dimensions by reducing it to a set of simple integro- 
differential equations. This work generalizes an approach we developed previously for hard discs. 
We numerically obtain both the pair correlation function and the virial coefficients for a fiuid of 
hyper-spheres in dimensions d = 4, 6 and 8, and find good agreement with available exact results 
and Monte-Carlo simulations. This paper confirms the alternating character of the virial series for 
d > 6, and provides the first evidence for an alternating character for d — 4. Moreover, we show 
that this sign alternation is due to the existence of a branch point on the negative real axis. It 
is this branch point that determines the radius of convergence of the virial series, whose value we 
determine explicitly for d — 4,6, 8. Our results complement, and are consistent with, a recent study 
in odd dimensions [R.D. Rohrmann et al., J. Chem. Phys. 129, 014510 (2008)]. 



I. INTRODUCTION 

In recent years, there has been much interest in the 
properties of hard sphere fluids. The study of such 
systems has played a central role in the understand- 
ing of classical fluids and serves as a starting point for 
the construction of perturbation theories of fluid prop- 
erties. A particular interest has been systems of hard 
hyperspheresii^i^'^'^'i'^'^iiSiiiiiiSiii^ (i.e. the generaliza- 
tion of spheres to dimensions larger than three, d > 3). 
There are several reasons for this interest. Firstly, as 
in a variety of interacting many-body systemSii^, one ex- 
pects studies of hard-sphere packings in high dimensions 
to yield great insight into the corresponding phenomena 
in lower dimensions, such as ground and glassy states of 
matter— li^. Secondly, it is hoped"'^'^ that analytical in- 
vestigations of hard hyperspheres in large spatial dimen- 
sions can serve as an organizing device for a systematic 
expansion in inverse powers of the dimension, d. From 
a completely different perspective, hard hypersphere sys- 
tems play an important role in communication theory. 
For example, it is known that the optimal way of send- 
ing digital signals over noisy channels correspond to the 
densest sphere packing in a high-dimensional space^i, so 
called "spherical codes" . 

It is common to describe hard-core fluids using approx- 
imate theories^^. There are two main reasons for using 
such approximate theories: on the one hand, a solution of 
the full problem was, and still is, extremely difficult. On 
the other hand, approximate methods have led to very 
good predictions in the low density phase. Among the 
most widely used approximations is the Percus-Yevick 
(PY) equation for d-dimensional hard spheres^*, which 
is exact to first order in the density of the fluid^^, p. 

A great deal of progress has been made towards un- 
derstanding the solutions of the PY equation in odd 
dimensions. In principle, the work of Baxter— and 
Leutheusser— reduces the PY equation to a set of non- 
linear algebraic equations of order 2^'^"'^)/^ for d > 3. 



However, these can only be solved analytically for d < 7 
so that results in higher odd dimensions must be found 
numerically. A more complete survey of the literature for 
odd dimensions may be found elsewhereii. 

In even dimensions the situation is much less de- 
veloped. In two dimensions an approximate numeri- 
cal solution of the PY equation was found by Lado^S. 
Leutheusser— was able to fit many of Lado's results using 
an ansatz for the direct correlation function. Other re- 
sults available in the literature for hard discs are solutions 
of the full problem based on Molecular Dynamics (MD) 
or Monte Carlo (MC) methoda^i2ai2ii22i^. Recently, the 
authors solved the PY equation for hard discs^, by de- 
veloping a method that reduces the problem to a set of 
integral equations that are solved numerically without 
major difficulties. 

In larger even dimensions there are some MD 
simulations^i^i^, MC simulationai, MC calculations^!^ 
and a few analytical results for the low order virial 
coefficients^ii. On the analytical front, Rosenfeld^l gen- 
eralized Leutheusser's Ansata^^ to higher dimensions and 
compared the results with the analytical results in three 
and five dimensions. However, to our knowledge, the PY 
equation has not been solved in any even dimension apart 
from d = 2. 

In this paper, we solve the PY equation for some even 
dimensions (d = 4,6 and 8). We do so by generalizing 
our previous work^, which was based on techniques bor- 
rowed from the resolution of crack problems^ and uses 
some results from Baxter's classical method^ii^. The 
main difference from the hard sphere case (and any odd 
dimension in general) is that the problem of finding the 
total correlation function and the direct correlation func- 
tion are coupled. This means that the present analysis 
necessarily yields both correlation functions and there- 
fore provides the equation of state. The advantage of the 
current method over previous approaches is that it pro- 
vides all the quantities of interest as power-series in the 
density. Thus, questions like existence of negative virial 



coefficients in four dimensions^ and more generally the 
radius of convergence of the scries in large dimensionsii 
can be tackled. 



II. THE FERGUS- YEVICK APPROXIMATION 

The pair correlation function g{r) is related to the 
direct correlation function c(r) through the Ornstein- 
Zernike equation bjJ^i^^ 



h{r) = c{r) + p h{r')c{\r-r'\)dr' , (1) 

Jo 

where p is the particle number density and 

hir) = 5(r) - 1 , (2) 

is the total correlation function'^^. The PY approxima- 
tion is a closure relation for Eq. ([1]) . For a hard-core pair 
interaction potential, this approximation reada^ 

g{r) = h{r) + 1 == 0, r < 1 (3) 

c{r) = 0, r > 1 (4) 

Interestingly, the PY approximation can be seen as a 
Random Phase Approximation (RPA) to some nonlinear 
field theory as shown previously^. 

Here and elsewhere, we take the diameter of the hyper- 
sphere to be unity. Thus in d dimensions (rf = 2(fc + 1), 
with fc > 0), we have 



which can be rewritten as 



1] 



^d(l/2) V^ 



fc+i 



(fc + l)!r/ 



where 
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(5) 



(6) 



is the volume of a d-dimensional hypersphere of radius 
R, rj is the packing fraction, and the space filling den- 
sity corresponds to 77 = 1. We define the d-dimensional 
Fourier Transform 

~f{q)^{2Trf+\-^ r'^+'Mqr)f{r)dr, (7) 

Jq 

where Jk is the Bessel function and the inverse Fourier 
transform for /(r) and Fk{r) are given by 

/•oo 

/(r) = (2^)-('=+l)r-M q'+'Jkiqr)fiq)dq, (8) 
Jo 

Applying the d-dimensional Fourier transform to 
Eq. 11]) yields 



h{q) = c{q) + ph{q)c{q) 



(9) 



[l~pc{Q)] l + ph{q) 



= 1. 



(10) 



Finally, the static structure factor s{q) of wavenumber 
q is related to the pair correlation function through 



s{q) = I + ph{q) = 



1 



1 - pc{q) 



(11) 



III. RESOLUTION 



The condition (j4]) together with fS]) imposes that c{q) 
can be written without loss of generality as 



c{q) 



fc+i 



t''+--J,+ i{qt)<t>{t)dt (12) 



where 0(t) is a real function. Substituting (fT2|) into ([8]) 
and simplifying by using the integral 6.575-1 in'^-, we 
find that c{r) can be expressed as 



c{r) 



^ t(j)(t) dt 

r 7P 



r2 IT 



<r< 1 



(13) 



Since c{r) is discontinuous yet finite at r = 1, one has 



m-ii-t") 



2N-1/2 



i-> r 



(14) 



In the following, we use the formulation of Baxter for 
the odd dimensional case, with the important difference 
that instead of solving for c(r) directly, we solve for (t>{t), 
and obtain c{r) using Eq. (fTS]) . We use the Wiener-Hopf 
method by defining 



Mq) 



<q) 



i-p 



27r 



fe+i 



t''+^Jk+iiqt)(t>{t)dt 

(15) 

One sees that A{q) — A{—q), A{q) ^ 1 as q ^ 00 and 
that A{q) has the same properties as the corresponding 
function defined in the odd dimensional case: it has nei- 
ther zeros nor poles on the real axis, since by definition 
s{q) has neither zeros nor poles for all g's. Therefore one 
can use the Wiener-Hopf decomposition of Baxter?ii22, 



A{q) = Q{q)Q{-q) 



(16) 



where Q{q) is analytic for 5(<7) > 0. Following the same 
steps as i n^ii28,35 ^j^g ^^.^^ show that Q{q) can be written 
as 



Q{q) = 1 - A / Q{t)e'^'dt 
Jo 

where A is a parameter defined by 
A = (27r)V 



(17) 



(18) 



and thus using ((T5l) - ((T7| 



27r 



s'^+^Jfe+i(9s)0(s)ds 



Jo Jo 

- \^ f ds f ds'Q{s)Q{s')e"i^'"''^ (19) 



Multiplying by exp(— i^/i), with < t < 1 and integrating 
with respect to q from — oo to oo gives 



'- / dss^+i4,{s) / dqq-^''+^\jk+^{qs) cos qt 
= Q{t) - \ I Q{s)Q{s - t)ds (20) 

which can be simplified to give 



{s'~t^Y(l){s)sds = 2''k\ Q{t)~X Q{s)Q{s~t)ds 

(21)' 

Differentiating fc-times with respect to t^ and once with 
respect to t gives 



m = (-1) 



k+l 



fc+1 



tdtj 



Q{t)-X Q{s)Q{s-t)ds 

(22) 
which is valid for < t < 1. Therefore, once Q{t) is 
known, (j}{t) is given by Eq. (|22p and c{r) is given by 
Eq. ([T3D. 

Now let us work on the function h{r). Since h{q) is an 
even function one can write without loss of generality 



K,) - ( - 
q 



1"+'^ Jf^+^{qt)il){t)dt (23) 



where ^(i) is a real function. Then, using Eq. ([8|) h{r) 
can be written in terms of ip{t) as 



h{T) 



tipit) dt 



r>0. 



(24) 



Combining this equation together with the condition ([3]), 



one obtains 

'■^ #(0 dt 



#(t) dt 

1 %/t2 - r2 TT 



0< r < 1. 

(25) 

This is an integral equation of Abel type that we can 
invert. As shown in'^^ the inversion of the equation is 
given by 



i,(t) = 



VT^t2 



— T^SV' S — 

S^ — r^ TT 



< i < 1 



(26) 



Eq. ([26|) is an integral equation that determines ?/'(i) for 
< t < 1 as function of t\j(t) for i > 1. Also, note that 
the behavior of 'i\)(t) near i = 1 is 



V(i)~(l-i2) 



2\-l/2 



as i — > 1" 



(27) 



Substituting the results of Eqs. HI]), HSll-dni) and ^ 
into Eq. (fTU)) gives 



Q(-9) 



1 - A / g(s)e*«"ds 





(28) 



1 + P — 
9 



t^+^Jk+^{qt)ip{t)dt 



Multiplying by exp(— igi) with i > and integrating with 
respect to q from — oo to oo we obtain 



2''k\ Q{s)S{s-t)ds- (s^ _ t^) V(s)sd(i29) 

+ A / dsQ(s) / (s'2-(s-i)2)fe^(s')s'ds' = 

JO "'It-s 

which leads to 

/oo 
(s^ - t2)feT/,(s)sds (30) 

nl /"OO 

= -A / dsQ(s) / (s'2 - (t - s)2) V(s')s'ds' , 

Jo J\t-s\ 

where 8(x) is the Heaviside function. Recall that ip{t) is 
defined only for i > 0. Therefore one has 



2'=fc! Q{t) - / (s2 - i2)fe^(s)srfs = _A / ds Q{s) / (s'^ - (t - s)^) V(s')s'f^s' < t < 1 (31) 

Jt Jo J\t~s\ 

D /•! /'OO 

{s^-t^)^i;{s)sds ^ X dsQ{s) {s'^-{t-sf)''i;{s')s'ds' t>l (32) 

't -/o J{t-s) 



r 



Replacing ■0(0 fo'' < t < 1 in Eq. ([3T|) with its value as given by Eq. ((26|) yields an equation that depends only 



on iplt) for i > 1, namely 



Q{t) = A{t) + / B(s, t)iP{s)sds - X 

Ji Jo 



A{t 



B{s',t~s)'iP{s')s'ds' 



Q{s)ds < t < 1 



(33) 



r 



with 



A{t) = 



B{s,t) 



(2fc + l)!! 



(l-t^y 



2Nfc+i 



2^'k\ 



i^' 



eti 



1 1 1 



i2'2'2 



(34) 
k ] (35) 



where /(z; a, 6) is the regularized Beta functional. On the 
other hand, using Eq. (j26p and differentiating fc-times 
with respect to t^ and once with respect to t one can 
rewrite Eq. ((32|l as 



V;(i) 



V;(i) 



-(-1)'M73T 



tdt 



-(-1)' 



fc+i ^1 



\ f d 



2'=fc! Vtdi 



,fe A /d 



"^ ^^ 2fcfc! Udi 



A{t-s)+ / B(s',t-s)V'(s')'s'^s' 



fc+i ^t-i 




k+l „1 



(s''-(t-s)")"^V(s')s'(is' 



(t-s) 







{s'^-{t-sYYij{s')s'ds' 



its) 



Q{s)ds 

Q{s)ds 1 < t < 2 (36) 

Q{s)ds t > 2 (37) 



r 



Our approach has reduced the PY problem for hard 
hyperspheres to the solution of the set of one-dimensional 
integro-differential equations ([55]) . ([5S)) and ([57)1 for the 
auxiliary functions ^(s) and Q{s). Once these functions 
have been determined, the physically relevant functions 
may be determined; g{r) from ([2]) and ([24]) . c{r) from 
P^ and P^ . We note that unlike the odd dimensional 
case-^, in even dimensions it is not possible to separate 
the problem of finding the direct correlation function c{r) 
from that of finding the pair correlation function g{r). 
This is because the behavior of the auxiliary function ■ip{t) 
for < t < 1 is coupled to its behavior for t > 1 through 
Eq. P^ . Although we were unable to find an analytical 
solution valid for all p, a numerical algorithm to find 
the numerical solution of these equations can easily be 
implemented. Before dealing with the numerical analysis, 
let us first consider the equation of state in the present 
formulation of the problem. 



IV. EQUATION OF STATE 

There are two methods used to calculate the equation 
of state when the radial distribution function, g{r), is 
known. Without the assumptions made in deriving the 
PY equationi^, these two methods would yield the same 
equation of state. The difference in the equations of state 
calculated using these two methods therefore provides an 



estimation of the error made by using the PY approxi- 
mation. The first method of calculating the equation 
of state uses the isothermal compressibility kt which is 
given by^ 

dp 



1 



p(3 ^KT = - 



siq = 0) , 



(3 \dP('='>yj. 

where (3 is the inverse temperature l/fcsT. 
Eqs. ([15]) and p2|) . it is easy to deduce that 

1 ^, 2(2^)V '•1 
s(0) (2fc+l)!! 



(38) 
Using 



i2('=+i)0(i)cit, (39) 



which can be further simplified by successive integration 
by parts to give 



1 



1 



2tt^P 
k\ 



dt 



{s'-t'Y^{s)sds 



(40) 



Using Eq. pT|) we may eliminate 0(s) to yield 



Q{t)-\ / Q{s)Q{s-t)ds 



dt (41) 



-!- = l-2A 
m JO 

Thus, the compressibility equation of state can be written 
as 



/3P(= 



- 1 



Q(i)-A' / Q{s)Q{s-t)ds 



(42) 
dt )■ X'dX' , 



The second method to obtain the equation of state is 
derived from the virial theorem and is given by 



/JPC) =p 



^k+l p2 



5(1+) 



(fc + 1)! 2 
FinaUy, using Eqs. Q and ^, Eq. (gSl) becomes 



1 



ttA 



where A is defined in (fTS]) . 



stl>{s) ds 

Vs^ - f TT 



(43) 



(44) 



NUMERICAL PROCEDURE 



As in our earher work2^, we solve Eqs. ([55|) . (|36p and 
P7p for the auxihary functions Q{s) and '^(■s) iteratively. 



However, the present method of solution is simpler than 
the one used previously, as explained below. We pose 
power series 

oo 

Q{t) = ^V9,(t), 0<t<l, (45) 

oo 

^{t) = ^AV,(i), t>\, (46) 

for the unknown functions, and substitute these power 
series into Eqs. ([55]) . (|55|) and (|57)) . At zeroth order, we 
obtain 



(7o(t) = v4(i); V'o(0 = 0- 
For i > 0, we find 



(47) 



'0*+i(i) 



-(-1)' 



tdtj 
(-1)'= / d 






i=o 



A(^ — s)Jj^i + / B{s' ,t ~ s)ipi-j{s')s'ds' 









(t-s) 



(s'"-(i-s)")>._j(s')s'ds 



r-l 



i=o' 



%+i(0 



/ B{s,t)^i+i{s)sds~y2 



j4(t — s)(5j^i + / B{s',t — s)tpi-j{s')s'ds' 



qj{s)ds 
qj{s)ds, l<i<2; (48) 

qj{s)ds, t>2; (49) 

gj(s)ds, < i< 1.(50) 



r 



We use an iterative procedure, starting with (|47|) . to and 
calculate successively ipi{t) and then qi{t) using (|i5|) - ([5(I)) . 
In the above formulation of the problem, the only diffi- 
culty is that one has to differentiate {k + l)-times. How- 



B. 



(c) 



B. 



(v) 



-fe+l 



2(A: + 1)! 



(53) 



ever, this is balanced by the fact that the integrals to be ^-^^^^ ^^^^^ coefficients must be found numerically. For 
computed have no singular behavior, in contrast with our 
earlier method of solution for the case A; = (i.e. d = 2)^. 
Rather than computing the equation of state for var- 



ious densities p (as in earlier work^^), we compute the 
virial coefhcients, Bi, which are defined by 



I > 1 we find 

.,) _ 2(2^)^('+i) 



B 



i+2 



i + 2 



pp^j2^^p'- 



(51) 



/ Qiit)-^ qj{s)qi^i^j{s-t)ds 

Jo „_n Jt 



J=0' 



Thanks to the iterative procedure presented above, these ^ 
coefficients, namely i3j , from the compressibility route 

(|^ . and i?j , from the virial route dU]), are directly 
given by the numerical resolution of the problem. 

The first two coefficients are identical and may be com- 
puted analytically: 



B 



(,) (27r)fe('+i) 



i+2 



2'=+i(fc + l)!7i V^^^T 



ds . 



(54) 
dt. 

(55) 
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(c) 



B 



(v) 



1, 



(52) 



The calculation of the first two virial coefficients via 
the two routes (c and v) produce identical results (coin- 
cident with the exact results) because the PY approxi- 
mation is exact up to first order— in p. What is not so 



9v/3 I 10 



27r 
27\/3 



Tt' (n _ 81^ I 38848 \ 

1728 \ IOtt "•" 15757r^ j 

r^^ /r, _ 2511\/3 I 17605024 

0592 Y' 2807r """ 606375^ 



2.062 

2.281 

0.576 

-0.021 



TABLE I: The exact expressions for B^ as well as their corre- 
sponding numerical values for some even dimensions obtained 
recently^. 



evident is that for the same reason both the third virial 
coefficient and the first order pair correlation function, 
g\ (r) , are also reproduced exactly by the PY theory. B3 
is given in closed form by^ 




B 



(c) 



B. 



(v) 



^2k+2 



Bl 



^3/4(^ + 1,5) 



B{k + %,\) 



2[(fc + l)!]2 

(56) 

where B{a, h) — T{a)T{b)lT{a + 6) is the beta function, 
Bx{a,b) is the incomplete beta function^"*. Note that 
for each integer dimension this expression can be written 
in a simpler form without the need for transcendental 
functions. However, it is not possible to write a simpler 
general form. Furthermore, gi{r) is given by 



9i{r) = l + 



tipi{t) dt 



VW^ 



= e(r-l)[l + 2%(r;l)] , 

(57) 
where a2(r; 1) is the scaled intersection volum o^^'^^ . 

It is also worthwhile mentioning that ,64 was recently^ 
evaluated exactly for all dimensions up to d = 12, though 
a closed-form formula for any d is not available. We 
reproduce the values of B^ for even dimensions d < 8, 
of interest here, in Table [D 

We note that the computation of the correlation func- 
tion g{r) and the equation of state do not present any 
significant difficulties. The integrands that must be com- 
puted in the calculation of the coefficients i?j , see ([55]). 
and g(r), see ([M]) . have weak square root singularities 
that may be dealt with by integration by parts or by 
subtraction of the singularity. We validate our numerical 
implementation of the iterative procedure outlined above 
by comparing our results with the earlier results for the 
case k = Q (i.e. d = 2)^^. 



VI. NUMERICAL RESULTS 



FIG. 1: (Color online) The correlation function g{r) com- 
puted from the PY equation (curves) and from Monte Carlo 
simulations^ (symbols) in d = 4. Results are plotted for 
p = 0.1 (A), p = 0.3 (O) and p = 0.5 (D). The PY re- 
sults are plotted taking 10 terms of the series (I45|l - (|46|l and 
were computed with Ar = 10"'^. 



A. Comparison with Monte-Carlo simulations 

We begin by comparing the pair correlation function, 
g{r), obtained from our solution of the PY equation with 
that obtained in the Monte-Carlo (MC) simulations of 
others^. For simplicity we present here numerical results 
obtained using the first ten terms of the series (|45p - (P5)) 
and a spatial resolution Ar — \Q~^. These are represen- 
tative of the results obtained by taking more terms and 
using a higher spatial resolution. 

Figures [2 [21 and [3] show g{r) for d = 4, 6 and 8, respec- 
tively with several different values of p. (The values of 
p used here are chosen to be smaller than the radius of 
convergence of the virial series, see ijVI Bl and WlCl be- 
low. For larger values of p, the series (|^5|) - P5)) does not 
converge, though it may perhaps be resummed to im- 
prove the convergence for larger densities.) We see that 
there is generally very good agreement in each case be- 
tween the PY results (solid curves) and MC simulations 
(symbols). We note that, as expected, the discrepancy 
is largest for r « 1 — similar to the behavior observed 
in other dimensions22i^. Interestingly, this discrepancy 
seems to diminish as the dimension increases. The good 
agreement between the MC results and our solutions of 
the PY equation shows that our method of solution works 
well. 



In this section, we present the results of our numerical 
computations in dimensions d = 4, 6,8. A brief descrip- 
tion of our numerical scheme is presented in Appendix 

El 



Virial coefficients 



We computed the virial coefficients using the two ex- 
pressions ([51)1 and ([55)1 for d = 2 as a check of our nu- 
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FIG. 2: (Color online) The correlation function g{r) com- 
puted from the PY equation (curves) and from Monte Carlo 
simulations^ (symbols) in d = 6. Results are plotted for 
p = 0.1 (A), p = 0.2 (0) and p = 0.3 (O)- The PY re- 
sults are plotted taking 10 terms of the series (I45|l - (|46|l and 
were computed with Ar = 10^^. 



1.15 



D) 




1.05 



FIG. 3: (Color online) The correlation function g{r) com- 
puted from the PY equation (curves) and from Monte Carlo 
simulations^ (symbols) in d = 8. Results are plotted for 
p = 0.1 (A), p = 0.2 (0) and p = 0.3 (O)- The PY re- 
sults are plotted taking 10 terms of the series (|45|l - (|46|l and 
were computed with Ar = 10~^. 



merical scheme. (Note that this is not a trivial verifica- 
tion since the numerical scheme is completely difi^erent 
to that used previously^^.) We obtain the same results 
as reported previouslyS^, though our estimation of errors 
leads to some slightly different values for the last digit of 
the higher order coefficients. 

The virial coefficients calculated by this method for 
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3.083 


3.083 


3.083 
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1.774 


2.466 


2.281 
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0.988 


1.602 


1.323 
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0.315 


0.875 


0.707 
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0.285 


0.463 


0.323 
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-0.103 


0.193 


0.161 
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0.278 


0.120 


0.061 


10 


-0.386 


0.0046 


0.038 


11 


0.642 


0.071 


— 


12 


-1.058 


-0.076 


— 


13 


1.793 


0.134 


— 


14 


-3.088 


-0.210 


— 


15 


5.402 


0.348 


— 



TABLE 11: Numerical values of the first fifteen virial coeffi- 
cients in four dimensions 



The B\ are the results from 



[v) 



Monte Carlo calculations presented previously^. Bl and 
Bl are the values found from the solution to the PY equation 
using Eqs. (|55p and (|54p respectively. 



B 



(v) 



B 



(c) 



B 



{MC) 



3 


2.276 


2.276 


2.276 


4 


0.189 


0.795 


0.576 


5 


0.51(7) 


0.342 


0.335 


6 


-0.72(2) 


-0.104 


-0.200 


7 


1.43(7) 


0.278 


0.389 


8 


-3.0(4) 


-0.509 


-0.688 


9 


6.8(6) 


1.043 


1.326 


10 


-16.(3) 


-2.255 


-2.696 


11 


40.(5) 


5.097 


— 


12 


-10(2) 


-11.94(7) 


— 


13 


26(5) 


28.86(1) 


— 


14 


-7(10) 


-71.51(5) 


— 


15 


19(15) 


181.09(1) 


— 



TABLE III: Numerical values of the first fifteen virial coef- 
ficients in six dimensions. The B^ are the results from 
Monte Carlo calculations presented previousl)*^. -B, and 
Bj are the values found from the solution to the PY equation 
using Eqs. (|55p and (|54p respectively. 



d = 4, 6 and 8 are reproduced in tables [TIllIVI These ta- 
bles show the virial coefficients resulting from both the 
virial (v) and compressibility (c) routes, as well as the 
virial coefficients obtained by earlier MC calculations^. 
Note that all the virial coefficients up to B4 are known 
exactly (see Eq. ([55)1 and Table U above), and agree with 

(c) 

the MC calculations. In general, we see that B^ ' is a bet- 
ter estimator of the true virial coefficient, B^ , than 

B^l 

We seem to see the same trends as Rohrmann et alr^ 
regarding the way in which the PY virial coefficients 
bound the true virial coefficient for d > 7 (though our 



B 



(v) 



B 



(c) 



B. 



(MC) 



3 


0.966 


0.966 


0.966 


4 


-0.13(1) 


0.065 


-0.021 


5 


0.21(3) 


0.074 


0.126 


6 


-0.30(2) 


-0.078 


-0.155 


7 


0.5(1) 


0.11(5) 


0.239 


8 


-0.9(0) 


-0.19(0) 


-0.406 


9 


1.(7) 


0.33(6) 


0.747 


10 


-3.(4) 


-0.6(4) 


-1.466 


11 


7.(5) 


1.(3) 


— 



TABLE IV: Numerical values of the first eleven virial coeffi- 
cients in eight dimensions. The B^ are the results from 
Monte Carlo calculations presented previouslj'^. B^ and 
Bj are the values found from the solution to the PY equation 
using Eqs. (|55p and (|54p respectively. 




results reveal that the transition in behavior they ob- 
serve occurs already at rf = 6, rather than d = 7). In 

(v) (c) 

particular, we see that B^ ' < Bi < B] for even i > 4 
and Sf ^ > B, > B'f^ for odd i > 5. 

Note also the intermediate behavior observed 
previouslyii for d = 5 is observed already with d = 4, 
though the details are a Uttle different: we find that 
Bi < Sf ^ < 5,^") for odd i > 7 and B, > fif ^ > sf > 
for even j > 8 — though one would need to see more 
real virial coefficients to say this with more confidence. 



FIG. 4: (Color online) Domb-Sykes plot for the virial coeffi- 
cients in six dimensions. The ratios are based on the virial 
coefficients for P'-"^ (A) and P'"' (Q)- The dashed lines show 
the relevant linear fits at large i and the points at the inter- 
section with the y-axis shows the extrapolated limits and the 
associated error bars. 



d 


„(") 

'/conv 


'/conv 


2 


1.008 ±0.002 


1.01 ±0.002 


4 


0.146 ±8 X 10"'^ 


0.150 ±0.003 


6 


0.024 ± 7 X 10"'' 


0.024 ± 6 X 10"'' 


8 


0.0055 ±8 X 10~* 


0.0051 ± 5 X 10"* 



C. Convergence of virial series 

A question of considerable interest is the radius of con- 
vergence of the virial series. This question is closely re- 
lated to the nature of the singularity closest to the origin, 
which was addressed recently in odd dimensionsii. An 
estimate of this radius of convergence may be made by 
using the Domb-Sykes plot^i^. This is an extension of 
the ratio test in which the ratio between successive terms, 
Bi/Bi^i is plotted as a function of 1/i. Such a plot of- 
ten reveals a linear behavior, the intercept of which then 
provides an estimate of the radius of convergence of the 
series, Pconv, via 



Pc 



lim 



B, 



B, 



(58) 



An example of such a plot is shown in Fig. |4] for the virial 
coefficients in six dimensions. 

We used the Domb-Sykes plot to determine the radius 
of convergence of the series^ for d = 2, 4, 6 and 8. The 
results are given in table |V] as the radius of convergence 
for the series in -q, ryconv, along with estimates of the 
error in each case (Eq. ([5]) may be used to convert ?7conv 
to Pconv)- These results are also plotted in Fig. [5] and 
combined with the results for odd dimensions d < 13 
obtained by the analysis of Rohrmann et al}^ . These 



TABLE V: The radii of convergence for the virial series, for 
both the virial (v) and the compressibility (c) routes, as a 
function of the dimension d. 



results seem to confirm the assertion that as d — *■ oo , 



Vc 



(59) 



which was conjectured by Frisch and Percusi^ for the 
full problem. The Domb-Sykes plots for d = 4, 6, 8 have 
negative intercepts with the vertical axis and positive 
slopes there (see, Fig. |4l for example). This shows that, 
in these cases, the singularity that limits the radius of 
convergence of the virial series is a branch point on the 
negative real axis^^ii^. In contrast, the Domb-Sykes plot 
for two dimensions suggests that the relevant singularity 
is a pole on the positive real axis. 



VII. DISCUSSION 

In this paper we generalized a semi-analytic method 
to solve the PY equation for hard discs^ to the case of 
hard hyperspheres in even dimensions. The essence of 
this approach is a reduction of the PY equation to a set 
of integro-differential equations for two auxiliary func- 
tions Q{s) and -0(5) as given by Eqs. (|55|) . ([5^ and ([57)1 . 
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FIG. 5; (Color online) The radius of convergence for the virial 
series, ryconv, as a function of dimension, d. The radius of 
convergence in odd dimensions (x) is taken from the work 
of Rohrmann et al^. The results in even dimensions are 
found from the relevant Domb-Sykes plots and show estimates 
based on P*'^' (A) and P'"' {£)). The dashed line shows the 
relationship jyconv = 2^'', which is believed to describe the 
large d behavior of 77conv 



of stated (for three-dimensional hard spheres) to higher 
dimensions^, are inferior to the PY theory since they 
are not able to reproduce such a branch point singularity. 
This characteristic is also missing from far more elabo- 
rated equations of stateiSiiS,. 

Another interesting prediction of the current work is 
that the exact virial coefficients Bi for d = 4 may change 
sign if they are computed for sufficiently large i. We 
observe that in PY theory with d — A the virial coeffi- 
cients determined by both routes are negative for even 
i > 12. This suggests that the exact coefficient may also 
become negative for sufficiently large i and will hopefully 
motivate the calculation of further coefficients using the 
methods of Clisby and McCoy^. More generally the ques- 
tion of negative values of the virial coefficients in various 
dimensions is pertinent^i^iii. 

As explained above, this method allows, in princi- 
ple, calculations to arbitrary precision, and it could be 
interesting to obtain more virial coefficients by doing 
so. Such progress might also allow for a comparison of 
the pair correlation function with Molecular Dynamics 
simulations^!^ at large densities. 

It could also be interesting to apply the approach de- 
veloped here to polydisperse mixtures2i^'— , to sticky 
hard spheres (i.e., hard spheres with an adhesive short 
range interaction)"*^ and to much higher dimensions. 



The correlation functions and the equation of state can 
be determined easily from these auxiliary functions. We 
suggest an efficient iterative numerical method to solve 
these equations and determine the auxiliary functions in 
d — 4,6 and 8. Using this method we are able to deter- 
mine the values of the virial coefficients within the PY 
approximation and compare them with the first ten virial 
coefficients for the full problem^. We also obtain results 
for the pair correlation function which compare well with 
the available MC simulations^. 

The principal advantage of this approach is that it pro- 
vides directly the virial series, and so it yields the equa- 
tion of state for all values of p at the same time provided 
that p < Pconv This is in contrast with other approaches 
where each value of p requires a separate calculation^^. 
An important consequence is that we can study the con- 
vergence of the series. We have shown that the virial 
series predicted by the PY theory for hypersphcrcs in 
even dimensions has a branch point on the negative real 
axis r] = — jyconv for all d > 4. This is observed in both 
the virial and compressibility routes to determining the 
equation of state, similar to what happens in PY in odd 
dimensions*^. The position of this singularity as a func- 
tion of the the dimension d is also consistent with the 
conjecture by Frisch and PercusiS for the full problem in 
large dimensions. The successful prediction of the singu- 
larity supports the idea that the PY theory approaches 
the exact problem as the dimensionality increases. An 
important conclusion from this discussion is that semi- 
phenomenological equations of state, such as the gener- 
alizations of the celebrated Carnahan-Starling equation 
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APPENDIX A: NUMERICAL SCHEME 

We solve the system of equations (|T7|) -([5n | by discretiz- 
ing in space using steps of size Ar. Integrals are calcu- 
lated using the trapezoidal rule, which is first order ac- 
curate. To calculate derivatives we developed two differ- 
ent schemes: one which used forward differencing (first- 
order accuracy) and the other using central differencing 
(second-order accuracy) . The results obtained with these 
two differentiation schemes are consistent with one an- 
other. Iterations proceed from i = using the values for 
ipo{t) and qo{t) from P7)) to determine V'i(^)- The func- 
tion tpi (t) may then be used with the discretized version 
of ([50)1 to determine qi{t). This process is then repeated 
N times, corresponding to determining the first N terms 
in the series expansions of -0(0 and q{t). The results 
presented in this paper typically have N — f 0. 

The virial coefficients were determined by numerical 
integration, again using the trapezoidal rule. To deter- 
mine the values of these coefficients more accurately, we 
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performed computations at different spatial resolutions, 
i.e. with different values of Ar. Plotting the behavior 
of these numerically determined coefficients as a func- 
tion of Ar we then extrapolated the observed trend to 
Ar — 0. The observed trend was linear for the scheme 
using differentiation by forward differencing, as expected 
since both discretizations involve errors of order Ar. The 
second scheme (using differentiation by central differenc- 
ing) is more complicated since the first order error of 
integration is mixed with a second-order error in differ- 
entiation. Typically with this scheme we find errors that 
scale like Ar'^". In Fig. [6] below we show an example of 
the convergence using these two schemes. 



FIG. 6; (Color online) An example of the convergence of the 
virial coefficients as a function of Ar for _Bg with d — 6 using 
the two numerical schemes described in the text: forward 
differencing (□) and central differencing (O)- For both data 
sets the dashed lines are the result of the fitting procedure, 
with a slope of 1 for forward differencing, and with a slope of 
3/2 for central differencing. Note that the plot is on a log-log 
scale. 



To determine the values presented in tables HIllIVI we 
use the values extrapolated to Ar — from both nu- 
merical schemes. An estimate of the error introduced in 
this extrapolation procedure was obtained using the 95% 
confidence interval for the value at Ar = 0. The results 
presented in tables HIllIVI are correct to 3 decimal places 
or as otherwise indicated. 
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